# Replication code for 
# Dom, R., & Roger, L. 'Debt or alive: Burundi's Fiscal Response to Economic Sanctions'.
# International Studies Quarterly, 2020.

rm(list = ls()) # "clear"

# Load packages. ----
# To install them all, uncomment:
# install.packages(c("tidyverse", "haven", "lubridate", "vars", "gridExtra", "xts", "Cairo", "strucchange", "stargazer"))
library(tidyverse) # A large bunch of nice functions (data wrangling, graphs, etc.)
library(haven) # To import stata files
library(lubridate) # To handle dates
library(xts) # To handle time series data
library(vars) # VAR analysis (SVAR, VECM etc.)
library(gridExtra)
library(Cairo) # For good looking pictures
library(strucchange) # For Bai-Perron test for break points
library(stargazer)
source("functions.R") # Load user-written functions

# Parameters ----

roll.average <- FALSE

# Load data ----

dta <- read_dta(file = "dataset_exp.dta") # Stata file saved after all transformations applied in explore.do
if(roll.average == TRUE) dta$Grants_Cash <- rollmean(dta$Grants_Cash,k = 3,fill = NA)
dta.sample <- dta[-1,] # Remove first observation, because it's missing for Ddebt_... (because of differencing)
if(roll.average == TRUE) dta.sample <- dta.sample[-nrow(dta.sample),] # Remove last if NA because of rolling average

dta.xts <- xts(dta.sample, order.by = dta.sample$date,frequency = 12) # Same dataset as proper time series object


endog <- dta.sample[,c("Grants_Cash", "Ddebt_int_total","exp_cur_tot","revenue")] # Endogenous variables
exog <- dta.sample[,c("cbDummy","cbDummy2","pierre")] # Exogenous variables (trend and intercept are included via VAR command)

# Plot data ----
ts.facetplot(data = endog, # Function from functions.R
             variables = c("Grants_Cash", "Ddebt_int_total","exp_cur_tot","revenue"), 
             index = dta.xts)

# VAR ----
# Lag-length selection
VARselect(y = endog,lag.max = 6,exogen = exog,season = 12,type = "both")

# All criteria indicate inclusion of only one lag, so we do that.
basic.var <- VAR(y = endog,p = 2,exogen = exog,season= 12,type = "both")

# Specification tests
arch.test(basic.var)
serial.test(basic.var,type = "PT.asymptotic")
serial.test(basic.var,type = "BG")
serial.test(basic.var,type = "ES")
stab <- stability(basic.var)
#plot(stab)

# Results
#summary(basic.var)

# Impulse response functions
basic.irf <- irf(basic.var,
                 impulse = c("Grants_Cash"),
                 response = c("Grants_Cash","Ddebt_int_total","exp_cur_tot","revenue"), # Must include impulse as well
                 cumulative = TRUE,
                 ortho = TRUE)

normal.irf <- irf.normalise(irf = basic.irf) # User-written function to express IRFs in unit changes instead of standard deviations

interest <- which(normal.irf$impulse == "Grants_Cash" & normal.irf$response == "Ddebt_int_total" & normal.irf$step == 11 |
                   normal.irf$impulse == "Grants_Cash" & normal.irf$response == "exp_cur_tot"  & normal.irf$step == 11   )
normal.irf[interest,] # Cumulated IRFs to report in appendix

p1 <- norm.irf.plot(irf = normal.irf, 
              imp = "Grants_Cash", 
              imp.label = "Cash grants",
              resp = "Grants_Cash",
              resp.label = "Cash grants",
              invert = TRUE)

p2 <- norm.irf.plot(irf = normal.irf, 
              imp = "Grants_Cash", 
              imp.label = "Cash grants",
              resp = "Ddebt_int_total",
              resp.label = "Domestic borrowing",
              invert = TRUE) 

p3 <- norm.irf.plot(irf = normal.irf, 
              imp = "Grants_Cash", 
              imp.label = "Cash grants",
              resp = "exp_cur_tot",
              resp.label = "Recurrent expenditure",
              invert = TRUE)

p4 <- norm.irf.plot(irf = normal.irf, 
              imp = "Grants_Cash", 
              imp.label = "Cash grants",
              resp = "revenue",
              resp.label = "Total revenue",
              invert = TRUE)

irf.combine <- grid.arrange(p1,p2,p3,p4, ncol = 2)
rm(p1,p2,p3,p4)

ggsave(filename = "IRFs.jpeg", 
       plot = irf.combine, device="jpeg", width = 12,height = 8)

# Key IRF numbers for appendix table
# NOTE: Lower and upper bounds may differ slighty from those reported in the paper
# depending on the system and software versions employed. This is because of the stochasticity 
# that is inherent to the bootstrapping procedure and minor computational differences between systems. 
# The observed differences are qualitatively negligible.

key <- which(normal.irf$impulse == "Grants_Cash" & 
               (normal.irf$step == 1 | normal.irf$step == 11))

normal.irf[key,]

# Intervention analysis

dta.xts$trendnum <- as.numeric(dta.xts$trend)

intervention <- lm(formula = dta.xts$Ddebt_int_total ~ dta.xts$pierre + dta.xts$trendnum)
summary(intervention)  # suggests a mean shift of 21936, vv significant.

intervention.exp <- lm(formula = dta.xts$exp_cur_tot ~ dta.xts$pierre + dta.xts$trendnum)
summary(intervention.exp)  # suggests a mean shift of 21936, vv significant.

stargazer(intervention, intervention.exp, title="Intervention Analysis",
          dep.var.labels=c("Borrowing","Spending"),
          covariate.labels=c("Intervention","Trend",
                            "Intercept"),
          omit.stat=c("LL","ser","f"),
          label = "tab:intervention", digits = 2)

# Bai-Perron test for structural change

bp <- breakpoints(dta.xts$Ddebt_int_total ~ 1 + dta.xts$trendnum)
summary(bp) # Allowing for m=1 breakpoint, it identifies March 2015 

bp.exp <- breakpoints(dta.xts$exp_cur_tot ~ 1 + dta.xts$trendnum)
summary(bp.exp) # Allowing for m=1 breakpoint, it identifies March 2015 


# Granger Causality

endog <- dta.sample[,c("Grants_Cash", "Ddebt_int_total")] # Endogenous variables
exog <- dta.sample[,c("cbDummy","cbDummy2","pierre")] # Exogenous variables (trend and intercept are included via VAR command)

gc.var <- VAR(y = endog,p = 2,season= 12,type = "both",exogen = exog)
causality(gc.var,cause = "Grants_Cash")
causality(gc.var,cause = "Ddebt_int_total")

# # SVAR (equivalent to what underlies orthogonalised IRFs, repeated to obtain A ["Omega"] matrix)

Omega = matrix(c(1,NA,NA,NA,
                 0,1,NA,NA,
                 0,0,1,NA,
                 0,0,0,1),
               nrow = 4)

SVAR <- SVAR(x = basic.var,Amat = Omega)
SVAR$A
SVAR$Ase

